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In this work, we report a novel technique in lattice QCD for studying the high momentum- 
transfer region of nucleon form factors. These calculations could give important theoretical input 
<«_^ I to experiments, such as those of JLab's 12-GeV program and studies of nucleon deformation. There 

Cn . is an extensive history of form-factor calculations on the lattice, primarily with ground states for 

^.^' both the initial and final state. However, determining form factors at large momentum transfer 

j^ ' {Q^) has been difficult due to large statistical and systematic errors in this regime. We study the 

nucleon form factors using three pion masses with both quenched and 2-|-l-flavor anisotropic lattice 
configurations with Q^ as large as 6 GeV^. These form factors are further processed to obtain 
|/-\ ' transverse charge and magnetization densities across 2-dimensional impact-parameter space. Our 

approach can be applied to isotropic lattices and lattices with smaller lattice spacing to calculate 
even larger-Q^ form factors. 



C^ 



PACS numbers: 13.40.Gp, 12.38.Gc, 14.20.Dh 



I) 

r-| ; I. INTRODUCTION 

The structure of hadrons is revealed by their interactions with various probes in scattering experiments. For electro- 

^ ' irragnetic interactions that do not change the particle content, we describe them in terms of the elastic electromagnetic 

On ', form factors, whose dependence on transfer momentum (Q^) provides information about hadronic structure at difFer- 

ON ' ent scales. Many experimental studies of these nucleon form factors have been conducted. Recently, a Jefferson Lab 

t_^ I experiment using both a polarized target and longitudinally polarized beam (so called double-polarization) revealed 

a non-trivial momentum dependence for the ratio G^^/G\j. This contradicted previous results using the Rosenbluth 

l/~J ' separation method, which suggested fJ-pG^^ / G^j^j ~ 1. The apparent contradiction highlights the possibility of a systcm- 

^^ , atic correction associated with two-photon exchange, affecting the Rosenbluth separation method more significantly 

^-•^ ■ than double-polarization. (For details and further references, see the recent review articles: Refs. [iHSl-) 

. . I Measuring the form factors at higher Q^ will help us to understand hadrons and challenge models based on quantum 

^ ■ chromodynamics (QCD). Future experimental facilities, such as the 12-GeV upgrade at Jefferson Lab, will provide 

precision data at large values of Q^. Experimentally, it is easier to study proton form factors at larger Q^, but 

rather challenging to study the neutron sector due to its neutral net charge and the lack of a free-neutron target. 
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d ' Current data for the neutron Ge form factor only reaches to « 4 GeV^ transfer momentum. However, the 12-GeV 
upgrade at Jefferson Lab will provide precision form factors to around 10-18 GeV^ for both the proton and neutron. 
Theoretically, perturbative QCD should converge better at larger Q^. However, it fails to describe recent BaBar results 
for 7*7 — ?► 7r° over a wide range of transfer momenta 4 GeV^ < Q^ < 40 GeV^[j, Q. This leaves an opportunity for 
nonperturbative QCD approaches to extend their region of applicability to these high-Q^ regions. Lattice QCD is a 
perfect candidate for the job. 

Studies in the nonperturbative regime of QCD theory have been difficult without resorting to model-dependent 
calculations or making approximations due to the strong coupling at long distances. However, by discretizing space- 
time into a four-dimensional lattice with a fixed lattice spacing and volume, we are able to compute the path integral 
(in tcrnrs of discretizcd versions of the QCD Lagrangian and operators) directly via numerical integration, providing 
first-principles calculations of the consequences of QCD. The techniques of lattice QCD have been applied to such 
varied phenomena as the spectroscopy of heavy-quark hadrons, the tower of excited baryon states, flavor physics 
involving the CKM matrix, hadron decay constants and baryon axial couplings. Using lattice QCD to study hadronic 
form factors will serve as valuable theoretical input for understanding hadronic structure within the QCD theory of 
the Standard Model. 

The nucleon form factors have been calculated on the lattice by many groups, and calculations are still ongoing [6|-19|. 
Recently, lattice calculations have also been used to calculate transition form factors involving excited nucleons 20[ . 



However, the typical Q^ range in lattice calculations of hadron form factors is less than 2.5 GeV . When one attempts 
higher-Q^ calculations, they suffer from poor signal-to-noise ratio; for example, see the case study for the pion in 
Ref. [m. 

In this work, we start to explore the possibility of reaching higher transfer momenta using lattice QCD by re- 
examining the conventional approach. To calculate the form factors, we need to calculate the three-point Green 
function, which requires an additional inversion of the fermion matrix over the whole lattice (with rank 12 x L'^ x T). 
Since traditionally this type of operation is resource intensive, especially for light quark masses, the parameters 
arc carefully tuned to maximize overlap with the ground-state nucleon on a specific gauge ensemble, and the same 
parameters are used throughout the whole calculation. This simplifies the three-point correlator analysis, since one 
only needs to consider one state in a given plateau region. However, when one increases the momentum of the nucleon 
state, the original fixed parameters are no longer optimal; thus, the signal dies out quickly after adding only a few 
units of discrete momentum on the lattice. We propose to keep multiple operators in the calculation with parameters 
tuned for both ground and excited states, but we also extend our analysis to account for these additional states. As 
a result, we can extract the best signals at each transfer momentum. Since our analysis explicitly treats multiple 
excited states, the ground state is safe from contamination. Further steps and details are addressed in Sec. [ill We have 



previousl y d emonstrated the idea in an exploratory quenched lattice calculation [23 and extended it to a dynamical 
ensemble [231 . In this work, we use multiple quark masses and improved statistics, giving a complete analysis for three 
pion-mass ensembles of both quenched and dynamical (N f = 2-1-1) configurations generated by the Hadron Spectrum 
Collaboration (HSC)t24J, i25|. Note that even in this study, we may suffer large systematic error due to the coarseness 
of the lattices (a^ « 0.100 and 0.123 fm, respectively). As finer lattices are generated by various lattice groups, one 
can apply the same approach to reach even higher Q^ and significantly reduce systematic discretization errors. 

Since the elastic form factors contain information about the spatial structure of the nucleon, we can convert our 
data into a description of its charge and magnetization densities. Due to the relativistic effects of the transferred 
momentum on the wavefunction of the nucleon, we cannot use a simple three-dimensional Fourier transformation 
without some recourse to models. Instead, we use the model-independent formulation of Ref. [2^ in terms of densities 
in a two-dimensional plane transverse to an infinite- momentum boost. 

In this work, we concentrate on the electromagnetic properties of nucleons. The structure of this paper is as follows: 
In Sec. [ni we provide details concerning the configurations for both quenched and dynamical lattices. The operators 
and two- and three-point analysis procedure are given, as well as various checks of the method. We detail how we 
extract the electromagnetic form factors from lattice calculations. In Sec. Illli we discuss the momentum dependence 
of the form factors. We also extract the electric-charge radii, magnetic radii and magnetic moments for the nucleon 
and compare with other TV/ = 2-1-1 calculations. We extrapolate the form factors to the physical pion mass and 
compare the quenched results to dynamical. We use our results to describe the spatial dependence of these quantities 
in a model-independent way as transverse charge and magnetization densities. Our conclusions and ideas for future 
improvements to the calculation are presented in Sec. IIVI 



II. METHODOLOGY AND SETUP 

In this work, we report on a calculation of nucleon form factors using anisotropic lattices, including both dynamical 
AT^ = 2 -I- 1 and quenched A^/ = 0. 

The 2-|-l-flavor anisotropic lattices used in this calculation were generated by the Hadron Spectrum Collaboration 
(HSC)[2J, [25[. These lattices use Symanzik- improved gauge action with tree- level tadpole- improved coefficients, 
yielding a leading discretization error at 0(a\, a^, g^a^)- In the fermion sector, they have anisotropic clover action[27[; 
the gauge links in the fermion action are 3-dimensionally stout- link smeared with smearing weight p = 0.14 and 
Up = 2 iterations. The renormalized gauge and fermion anisotropics are around £, = 3.5 (that is, as = 3.5at), and 
the inverse of the spatial lattice spacing is about 1.6 GeV. For more details concerning the lattices and their action 
parameters, please see Ref. [2^. From the HSC ensembles, we use the 16"^ x 128 lattices with pion masses of 875, 
580 and 450 MeV. Quark propagators on the lattices are evaluated for source and sink operators with five Gaussian 
smearing parameters: a G {0.5,1.5,2.5,3.5,4.5}. Six, four and two time sources (respectively, from light to heavy 
pion mass) are used, and a total of around 200 configurations are used from each ensemble. The quark propagators 
are calculated under antiperiodic boundary conditions in the time direction, while the spatial ones remain periodic. 
We construct hadronic two-point correlators from all possible source-sink smearing combinations; however, for three- 
point correlators, we reduce the computational burden by keeping only the diagonal source-sink smearing-operator 
combinations. 

We also use quenched 16'^ x 64 lattices with anisotropy £ = 3, using Wilson gauge action with fi = 6.1 and 
stout-link smeared |28[| Sheikholeslami-Wohlert (SW) fermions|29j with smearing parameters {p^Up} = {0.22,2}. The 
parameter v is nonperturbatively tuned using the meson dispersion relation, and the clover coefficients are set to 




FIG. 1: Nucleon efTective-mass plots with (|p| =0 and 5 in units of 



from the m,^ — 450 MeV ensemble. The square, 



diamond and triangle points correspond to smeared-smeared correlators with a = 0.5, 1.5 and 4.5 respectively. 



their tadpole-improved values. The inverse spatial lattice spacing is about 2 GeV, as determined by the static- 
quark potential, and the simulated pion masses are about 480, 720 and 1100 MeV. In total, we use 400, 200 and 
200 configurations respectively at each pion mass. On the quenched lattices, we use only three Gaussian smearing 
parameters: a € {0.5,2.5,4.5}. For both two-point and three-point hadronic correlators, we calculate all 9 possible 
source-sink smearing combinations. 

We construct correlators with the quantum numbers of the nucleon using baryonic interpolating operators of the 
form 



X'^ix) 



^abc inT 



[qVix)C-f,q',{xMix\ 



(1) 



where C is the charge conjugation matrix, and qi and 52 are one of the quarks {u, d}. For example, in the case of the 
proton, we want qi = u and (72 = d. Two-point correlators are derived from these interpolating fields as 



where p is the baryon momentum, the spin projection F = "2 ^'^'^ ^ ^'-^"^ ^ index over the different smearing 
parameters. Eq. [2]can be decomposed in terms of energy eigenstates: 



T'IUt:p) = T. 



En + mn 

2E„ 



Zn.AZn.BS 



-E„(p)t 



(2) 



where n indexes over the basis of nucleon energy eigenstates. (These states are defined to be normalized as 
(0|(x^)^|p, s) = Zun{p, s) with nucleon spin-1/2 interpolating field x^)i the spinors in Euclidean space satisfy 



y^^UN{p,s)uN{p,s) 



E{p)j* — i-f • p + m 



(3) 



Conventionally, one would choose a smearing parameter that optimizes the signal of the zero-momentum ground- 
state two-point correlators and carry out the form-factor calculation at various transfer momenta. We present an 
example from the dynamical 450-MeV ensembles. The left-hand side of Fig. [l]shows the effective mass plot of the zero- 
momentum two-point nucleon correlator with all diagonal (with A = B) Gaussian smearing parameters. Examining 
the behavior of these correlators, one usually chooses the smearing parameter that contains the least excited-state 
signal, since this makes the analysis of the ground state simpler. In this example, the correlator with Gaussian 
smearing parameter of cr = 4.5 is a good candidate: the excited-state signals die out around t = 10, giving enough 
data points to extract ground-state form factors from t S [10, 30], if we put the sink around t = 40 (~ 1 fm source-sink 
separation). If we increase the magnitude of the momentum, say to |p|^ = 5 in units of -^a^^, we can see the signal 
from the a — 4.5 correlator decays significantly. This is not surprising, since this smearing parameter was chosen to 
"filter out" higher-energy contributions. So when increasing the momenta, the signal for the ground state starts to 
disappear from broadly smeared sources, while it remains clear for smaller values of a. To improve the quality of the 
signal at higher momenta, we should use multiple a and explicitly subtract any excited-state contributions (ideally 
more than one) to make sure that the ground state will be free from them. 

To extract form factors from three-point correlators, we need £"„ and Z^ as inputs from analyzing two-point 
correlators. Here, we apply the variational method[30| to extract the principal correlators corresponding to pure 
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FIG. 2: Summary of nucleon masses used by all currently existing Nf =2 + 1 nucleon electromagnetic form-factor calcula- 
tions fiiCi 



energy cigenstatcs from our matrix of correlators. The N x N Gaussian smeared-smeared correlation matrix (with 
N — 5 for dynamical and 3 for quenched) can be approximated as 



Ci. 



n=l 



tEn 



with eigenvalues 



An{t-, tr) 



-{t-tr.)E„ 



by solving the generalized eigensystem problem 



C{t)V = X{t,tr)C{tr)V. 



(4) 



(5) 



(6) 



where V is the matrix of eigenvectors and tr is a reference time slice. The resulting 5 eigenvalues (principal correlators) 
A„(i, tr) are then further analyzed to extract the energy levels E^. In practice, we are only interested in the lowest 
two; the extra higher states provide buffers against contamination by higher excited states due to the lack of the 
orthogonality in our operators. Since they have been projected onto pure cigenstatcs of the Hamiltonian, the principal 
correlator should be fit well by a single exponential and double checked for the consistency of the obtained energies. 
The leading contamination due to higher-lying states is another exponential having higher energy; we use a two-state 
fit to help remove this contamination. The overlap factors (Zn) between the interpolating operators and the n state 
are derived from the eigenvectors obtained in the variational method. Since we diagonalize the correlator matrix 
independently on every time slice, Z„ is a function of time; even though the time-dependence is mild, we choose the 
best value of Z„ by minimizing the difference between the correlator matrix reconstructed from Z and E and the 
original two-point correlator data. 

The nucleon masses from the dynamical ensembles used in this work are summarized in Fig. [2l also shown are other 
Nf = 2-1-1 nucleon masses used in published nucleon form-factor calculations. Figure [5] summarizes the pion and 
nucleon masses used by various groups who have calculated nucleon isovector Dirac and Pauli radii using Nf = 2 + 1. 
"AnisoClover" uses anisotropic lattice with clover actions; the pion masses ranges from 450-850 MeV with spatial 
lattice spacing and size around 0.123 and 2 fm. "RBC/DWF" carried out calculations using domain-wall fermions 
(DWF) on ensembles with pion masses of 330-670 MeV, spatial lattice spacing 0.114 fm and box size 2.7 fm|17|. 
"LHPC/DWF" also used DWF ensembles with the same spatial volume but at a smaller lattice spacing (0.084 fm) 
than "RBC/DWF", focusing on the pion-mass region 300-400 MeV[la|. "LHPC/Mixed" used a staggered- fermion 
sea, DWF valence with pion masses 290-760 MeV and a = 0.124 fm, L = 2.5 fm; they include one additional point 
with lattice size 3.5 fm for a 350-MeV pion[l^- In this work, we emphasize calculation of larger-Q^ quantities, 
necessitating the use of rather heavy quark mass inputs. In future calculations, we plan to have lighter quark masses 
and larger volumes. We can see from the plots that our nucleon masses fall nicely onto the trend outlined by other 
Nf = 2 + 1 ensembles. Due to the higher quark masses used in this calculation, we do not find noticeable finite- volume 
effects despite our small « 2 f m box size. 

To calculate the nucleon electromagnetic form factors, we first calculate the matrix element {x^ {Pf)\^'^\x^ (Pi)) > 
where V^ = q^f^q is the vector current with q being either an up or down quark, and P{ij} are the initial and final 



nucleon momenta. Wc integrate out the spatial dependence and project the baryonic spin, leaving a time-dependent 
three-point correlator of the form 



n n' 
X ^TapUn'{pf , s') p{Nn'{pf , s') \V^,\ Nn{pi, s))Un{Pt, s)a, 



(7) 



where fn.n'iPfiPi, E'^, En, t, ti,tf) contains kinematic factors involving the energy En and overlap factors (Zn) obtained 
in the two-point variational method, n and n' are the indices of different energy states and Zy is the vector-current 
renormalization constant (which is set to its nonperturbative value). The projection T used on the quenched lattices 
is Tmix ~ 5(1 + 74)(1 + 17573), and the dynamical lattices use T4 = |(1 -K 74) and T53 = 1(1 + 74) (17573) ■ The source- 
sink separation (tf — ti) is 34 and 39 time slices in lattice units on the quenched and dynamical lattices respectively. 
(These give source-sink separations about 1.13 and 1.36 fm.) We use four final momenta (p)' = ^{0,0, 0}a~^, 
^{— 1,0, 0}a~\ ^{-1, -l,0}a~^, ^{-2,0,0}a~^) and vary the initial momentum over ah p*; = ^{n^;, nj^,?T,2}a~^ 



with integer rix^y^z and n^ 



'Tiy + n', 



< 10. By fitting the time dependence of the three-point correlators to the form of 



Eq. [7]with n and n' restricted to and 1, we extract the ground state and matrix elements involving the first excited 
state. In this work, we only concentrate on the ground-state matrix element with n — n' = 0. 

With smeared fermion actions, it has been seen on three-flavor anisotropic lattices with tree-level tadpole- improved 
fermion-action coefficients that the nonperturbative coefficient conditions are automatically satisfied |24l |32|. Similar 
behavior has also been observed in another quenched study [33|, where the nonperturbative coefficients or renormal- 
ization constants in a smeared fermion action differed from tree-level values by a few percent. So the local vector 
current used here is 0{a) on-shell improved with the improved coefficient set to its tree-level value. 

The form factors are then extracted from the vector-current matrix elements for any nucleon state N through 



{N\V^\N)^{q) = UNip-f) 



Fliq^h, 



F2{q^) 
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(8) 



with q = pf — Pi. By imposing the same projection matrix T, vector-current matrix elements, (N \Vp,\ N) (with 
n = n' = 0) and momenta used in the lattice calculation in Eq. [3 we obtain a series of equations corresponding to 
the various lattice correlators. The overdetermined system of linear equations allows solution for the Dirac and Pauli 
form factors i^i,2- 

To make sure our analysis is correct, we check our procedures against the traditional one, the "ratio" method, 
where one takes ratios of three- and two-point correlators: 
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(9) 



with P and G corresponding to a point source and a Gaussian-smeared source. (We can also use another Gaussian 
smearing to replace the point sources.) In this ratio, the exponential time dependence is canceled and one only needs 
to fit a constant. Note that this method only assumes that the region where one extracts the matrix elements contains 
only ground-state signal. Otherwise, the matrix element will contain the contribution of unwanted excited-nucleon 
states. Here we demonstrate the method using the data from the quenched calculation with 720-MeV pion mass, 
p)' = {0,0,0}. We use the largest Gaussian smearing {a — 4.5) three- and two-point data, which has good overlap 
with ground state and not very much excited-state signal in the two-point effective-mass. Fig. [3] shows that both 
methods have consistent d-quark contributions to the Dirac form factor at pion mass 720 MeV. The star points are 
obtained using the conventional ratio method, and the circular points arc from our analysis combining 3x3 three- 
point correlators. Note that we not only get consistent numbers for Fi, but at large momenta, our fitting approach 
dramatically improves the signal. This is because when one tries to improve the ground-state signal (normally done 
by examining the hadron effective-mass time dependence at zero momentum) with a choice of smearing operator, one 
wipes out not only the excited states but also the higher-momentum states. Therefore, when one tries to project 
onto these higher momenta, there is no doubt that the signal-to-noise ratio will worsen. Since we are explicitly 
considering excited states in our analysis, our sources which have good overlap with higher-momenta will be only 
lightly contaminated by excited-state signals. 
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FIG. 3: Unrenormalized d-quark contribution to the nucleon Dirac form factor Ff obtained from the ratio approach (stars) 
and the method adopted by this work (circles) on a quenched ensemble with m,r ~ 720 MeV 



III. NUMERICAL RESULTS 



In this section, wc compare the quenched and dynamical nucleon form factors at Q^ as large as 6 GcV^ and examine 
their dependence on the pion mass. We also consider the charge radii and magnetic moments; these quantities are 
commonly derived from the electromagnetic form factors near or at Q^ = 0. Using calculations on ensembles at 
different pion masses, we extrapolate the form factors and their ratios to the physical pion mass. Later in the section, 
we study the large-momentum region and look at the charge distribution. 

In Sec. ini we described how the Dirac (Fi) and Pauli (F2) form factors are obtained from lattice calculations. 
Another common set of form factor definitions, widely used in experiments, are the Sachs form factors; these can be 
related to the Dirac and Pauli form factors through 



Gm(Q') 






(10) 
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In this work, we only calculate the "connected" diagram, which means the inserted quark current is contracted 
with the valence quarks in the baryon interpolating fields. The "disconnected" contributions are notoriously noisy and 
difficult to calculate. Previous works have attempted to estimate their contribution through an indirect approach by 
applying "charge symmetry" to the connected ones (including the the octet family) [3J, Is^] ■ For G|^ (the strange con- 
tribution to the magnetic form factor) at small transfer momenta, Adelaide- JLab Collaboration got —0.046(19) using 
quenched lattice data with chiral perturbatively corrected vacuum [36| and another dynamical calculation obtained 
G%i = —0.06(3) using mixed action[33l with pion mass as low as 360 MeV. The same approach can be applied to the 
electric form factor, yielding G|;(Q'^ = 0.1 GeV^) = —0.0049(56) in the same dynamical calculation [33 ■ Recently, 
another direct calculation from xQCD Collaborationist using the Z4 stochastic method (with unbiased subtractions) 
on Nf = 2 + 1 clover fermion lattices, found G|^ = -0.015(23) and G%{Q^ = 0.1 GeV^) = 0.0022(19). 

Another independent work done by the Boston group using multigrid techniques on 2-flavor Wilson anisotropic 
lattices with 400 MeV pions also observed G|; consistent with zero over Q^ ranging 0.1-1 GeV^jS^- Both direct and 
indirect approaches indicate the disconnected contributions to the electromagnetic form factors are small, O(10~^), at 
low momentum transfer. One would expect the disconnected contributions to become yet smaller at large momentum 
transfer, where the strong coupling between the loop and the baryon becomes weaker. The light-quark disconnected 
contribution turns out to be only moderately larger than the strange one. Using the isospin symmetry approach, 
the ratio of the strange-quark disconnected contribution to the light-quark was estimated to be 0.16(4) and 0.14(4) 
for electric and magnetic form factors, respectively [Sa, \^. We are not aware of any direct calculations of the 
disconnected quark ratios in the electromagnetic form factors; however, they exist in many other quantities. For 
example, the strange-to-light ratio of the first moment of the quark momentum fraction is 0.88(7) [4l|, and the ratios 
for angular and orbital momentum contributions to the nucleon spin are also around 1[43. It should be safe to 
estimate that the light disconnected contribution to the electromagnetic form factors is on the same order as the 
strange disconnected ones: 0(10^^) at low momentum transfer. 



We expect the disconnected contributions to become yet smaller as the momentum transfer increases. From the 
mesonic cloud point of view, the disconnected diagram is suppressed by an additional factor of 1/Q^ relative to the 
connected diagrams; thus, its relative contribution is suppressed at large momentum transfer. For a form factor 
whose typical magnitude in the low- momentum region is 0(1) (such as proton form factors and neutron magnetic 
form factor), the disconnected contributions are smaller than the statistical errors at most momentum points and 
the systematics from sources such as the nonzero lattice spacing. However, the disconnected contributions could 
significantly affect small-magnitude form factors such as the neutron electric form factor. We nevertheless include 
results from this channel for comparison to other lattice results but remind readers that the neutron electric form 
factor has larger uncertainty than the other form factors in this calculation. 

Figures[3]and[5]summarize the quenched and dynamical calculations of the electric G^", magnetic G^}" and isovector 
Dirac and Pauli F^j form factors as functions of momentum transfer Q^; the different symbols (colors) indicate 
different quark masses (or equivalently, pion masses) used as input to the calculation, as indicated in the captions of 
the figures. The dashed lines on these plots are taken from a parametrization of the experimental data[l|, |43|. In the 
case of proton, the parametrization used should be valid up to around 6 GeV^ for Ge and 25 GeV^ for Gm- The 
neutron form factors are only accurately known to 4.5 GeV^ for G\[ and 1.5 GeV for G'%] we extrapolate all the 
experimental parametrizations through 6 GcV'^ for comparison. Later in this section, we will discuss the extrapolation 
of the lattice data to the physical pion mass and the form-factor ratios. 

From the experimental parametrization of the Sachs form factors G^^\j, we can obtain Dirac and Pauli form factors 
-F'1,2 by reversing the definitions in Eq. 1101 

F^m = ^£±1^ (12) 

1 + r 
F.iQ') = ^%±^ (13) 

l+T 

for both proton and neutron form factors, where t = Q'^/(2A/)^. The isovector F"" form factors are just F'p — F". We 
can further estimate the individual quark contributions F^^ i^i terms of Ff'2 as 

FUQ') = Ph + '^Fla (14) 

FtAQ^) = 2FI2 + FI2. (15) 

All the form factors are nonperturbatively renormalized such that G'^(O) = 1. The largest available transfer 
momentum from the dynamical lattice is limited by the coarse lattice spacing. The inverse of the spatial lattice 
spacing is only around 1.6 GeV, while in the quenched case it is 2 GeV. Application of the same calculation procedure 
on finer lattices could easily extend the largest available transfer momentum. The pion masses used are 1080, 720 and 
480 MeV for the quenched ensembles and 875, 580 and 450 MeV for the dynamical ensembles. Both quenched and 
dynamical form factors display a trend toward the experimental parametrization as the pion mass decreases except 
for G^. The quenched form factors have milder pion-mass dependence with around 300-MeV separation, compared 
with the dynamical calculation. The lightest pion masses in the quenched and dynamical calculations are roughly 
the same, allowing us to see how the sea-quark contribution infiuences the form factors. It is immediately obvious in 
G|, that the addition of the sea-quark degrees of freedom lowers the form factors toward the experimental line. G^ 
has the largest noise-to-signal; the magnitude is at the order of the disconnected contributions mentioned earlier this 
subsection which could be influencing the connected lattice points significantly. 

The lattice data show some broad trends: As the lattice momentum used in any given point increases, the noise 
also increases. In this calculation, this means that the highest transfer-momentum points will have larger uncertainty, 
but also that there will be some points with smaller transfer momentum that have large uncertainty simply because 
they are constructed using large initial and final momenta. For purposes of clarity, we omit some of these points with 
very large errors; all points are included in the fits, although points with large error have little infiuence. 

A. Isovector Radii 

The size of the nucleon characterized by the effective charge and magnetic radii can be determined from the 
electromagnetic form factors. We first examine commonly calculated quantities using Dirac and Pauli isovector 
currents, where the disconnected contribution is highly suppressed due to isospin symmetry. 

The isovector Dirac and Pauli mean-squared radii can be extracted from the isovector electric form factors F'12 via 
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FIG. 4: Nucleon form factors on quenched configurations with pion masses of 480 (brown), 720 (orange) and 1080 (purple) 
MeV. The dashed fines are a parametrization of experimental form-factor data[l|, |43] 



Most groups have studied radii with the Q^ dependence over ranges 0.5-2.0 GeV^ among their own data and found 
the extracted radii to be independent (within the statistical error bars) of Q^ choice [iSl [la. [l8| . 

It is commonly agreed that a dipole extrapolation should be used for F^ , but opinions differ as to whether a dipole 
(a((5^ + &)~^) or tripole (a(Q^ + 6)"^) is preferred for F^. Refs. [ia[l3 found insignificant differences between results 
for either choice, while Ref. [1^ observed some discrepancy and adopted the numbers from the tripole for (r|). A 
summary of all the Nf = 2 + 1 lattice calculations of the isovector radii can be found in Fig. [HI Note that only the 
statistical errors are shown in this figure. 

The Dirac and Pauli mean-squared radii from the dynamical ensembles are summarized in Fig. |6] along with other 
Nf = 2 + 1 lattice calculations and lowest-order heavy-baryon chiral perturbation theory (HBXPT) using experimental 
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FIG. 5: Nucleon form factors on 2+1 configurations with pion masses of 450 (brown), 580 (orange) and 875 (purple) MeV. The 
dashed lines are a parametrization of experimental form- factor data[3, 143 ] 

inputs [4J|. Our results are nicely in agreement with isotropic TV/ = 2+1 calculations having various sea and valence 
fcrmion actions; this demonstrates the universality of the lattice QCD calculations. 



B. Magnetic Moments 



Calculations in a finite volume (without using twisted boundary conditions or an external magnetic field) cannot 
jive us the Q^ « Pauli F2 or magnetic Gm form factor. However, we are interested in the (anomalous) magnetic 
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FIG. 6: Summary of the isovector Dirac and Pauli mean-squared radii from all currently existing Nf = 2 + 1 nucleon electro- 
magnetic form- factor calculations [l5l-[la| . The dashed line indicates the leading-order HBXPT prediction. 



moments: 
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0) 
= 0), 



(17) 
(18) 



where x can be the proton, neutron, or up or down quark flavors. Two approaches are considered: The first method 
is to perform a dipole fit to the Q^-dependence of the form factors. However, due to the limited minimum momentum 
available (a constraint related to the lattice box size), the fit form can be poorly constrained in the near-zero Q^ 
region, resulting in fit-instability in obtaining the magnetic moments. The second method can be applied in the small 
Q^ region by polynomial fitting to the ratio of the magnetic to the electric form factor, Gm/Ge (or F2/F1). The two 
form factors are expected to be functions of r = 3^-3-- When Q^ is small, we can Taylor expand both form factors 
(except for the neutron): 



Gm/Ge = a^/Ge(0)(1 + air) + 0(t2) 
F2/F1 = k/Fi(0)(1 + air) + 0(r2). 



(19) 
(20) 



Both methods (dipole fit and linear fit to form-factor ratios) have been performed in a previous study on most of the 
octet baryons[18| and the latter method was found to have good consistency and improved stability. The magnetic 
moments are related via Hx = Kx + Bx, where e is the electric charge and x can be proton or up or down quark flavors. 
Here we will only show the dipole form to extract the anomalous magnetic moments from isovector form factors; 
both approaches give results consistent within statistical errors. To better compare our results with other studies, 

we convert k into the natural units of the nuclear magneton (/xjv — 2;^) by multiplying by — rs^. Our dynamical 

results are displayed in Fig. [7] along with other Nf = 2 + 1 results. The k obtained from tripole fitting are consistent 



with dipole results; similar observations were made by other Nf = 2 + 1 calculations [15l4l7| . Once again, we observe 
the universal behavior among different groups and lattice-parameter choices. We also note that the lattice results 
display a rather mild dependence on the pion mass in the range from 300 MeV to 850 MeV, yet they are only 2/3 
the experimental value. We may again expect a rapid rise in the smaller-pion mass region if lattice QCD correctly 
reproduces the experimental values. In Ref. [18|, the magnetic moments of the octet baryons were studied with 
SU(3) next-to-leading-order (NLO) HBXPT formulae; they found large discrepancies from the lattice points (with 
the lightest pion at 350 McV), suggesting that NNLO effects are significant. We might have similar expectations for 
the isovector anomalous magnetic moments. 
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electromagnetic form- factor calculations [15l-[la| . including this work 



1 nucleon 



C. Chiral Extrapolation of Form Factors 

Phenomenologically, as done in Refs. [ll, |43, it is common to describe the Q^-dependence of the experimental Ge.m 
form-factor data using a dimcnsionlcss parameter r = .^^ via an expansion form 



where fc is a selected integer and Go is 
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(21) 



Note that in the continuum, there is no difference between 



extrapolating the data in terms of Q or t, since the tojv is fixed. However, in a lattice calculation where the sea 
and valence quark masses vary, m}^ is a function of the pion mass. Using a dimensionless parameter should help to 



-^uA 



smooth the chiral extrapolation. We could also extrapolate F^ ^ through a similar procedure as G^ ^j, except that 
F2 should have one fewer power of r; that is, 



F, = 



ao 



T,L? a^T' 



Fo 



1 + Eti&»^* 



(22) 



To extrapolate to the physical pion mass, we can use one of several different approaches: First, we could use a 
simultaneous fit to the r and TOtt dependence where the fit parameters, a^ and bi in Eqs. [21] and [22] are in terms of 
rn^. A second approach is to fit data for each lattice momentum as a function of Q^{m^). Each data point f{Q^, m^) 
composes with the same kinematic momentum combination can then be extrapolated to the physical point, where we 
can further apply Eqs. [2l1 and [22] to fit the Q^ dependence. 

In this work, we simultaneously fit the ra-^ and Q^ (or r) dependence of the lattice data, expanding each fit parameter 
in terms of the pion mass: ai = al + a\ "ij. The pion masses on these ensembles are heavy enough that the linear 
ansatz for squared-pion-mass dependence should hold, simplifying the complicated extrapolation. This doubles the 
number of fitting parameters, but increases the number of data points by a factor of the number of pion masses used 
(3 in our case). Figure [5] shows an example of the extrapolations on the isovector Dirac and Pauli form factor in the 
dynamical ensembles. The x^/dof for F^2 ^'^e 1.4 and 2.0, respectively. The points and their corresponding symbols 
are the same as those in Fig. [5] and now we have the Q^ extrapolation line with statistical error band going through 
these lattice points. Further, the lowest line/band represents the extrapolated form factors at the physical pion mass. 
We repeat a similar process for both quenched {Nf — 0) and dynamical {Nj =2 + 1) on other form factors. 

We find good agreement for the extrapolation to the physical pion mass in Fig. [5] that is, the curvature (which 
corresponds to the charge radii) as a function of Q^ fits nicely with interpolating forms derived from experimental 
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FIG. 8: Isovector Dirac and Pauli form-factor extrapolation using Nf =2 + 1 calculations in this work. The shaded bands 
show the simultaneous fit at the pion mass corresponding to each data set, where the color of the band matches the color of 
the points. The lowest gray band is the extrapolation to the physical pion mass. 



data. This seems to contradict what we obtained in Subsec. IIII Al Note that the difference in terms of Q^ and r 
dependence is a factor of 4(77i^')^. This makes no difference in the continuum, since the nucleon mass is a fixed 
constant. However, in lattice-QCD calculations, where the quark masses are varied (as shown Fig. [2]), this could 
create a strong dependence. The good agreement in the radii extrapolation might be due to cancellation between the 
chiral curvatures in the nucleon mass and the radii. We compare our results with RBC/UKQCD Nf = 2+1 results 
and observe similar behavior. Figure |9] shows the same Nf = 2 + 1 products of the Dirac and Pauli radii with 4(to^')^ 
from our and RBC/UKQCD Nf = 2 + 1 results. The results are encouraging; the notorious charge radii problem may 
be ameliorated by looking at dimensionless quantities. In the Dirac radii product case, the discrepancy is improved 
to within a few standard deviations and even smaller for the Pauli radii product. However, further investigation and 
improvement in the statistical errors of these calculations will be required to fully understand the nature of these 
dimensionless quantities. 
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FIG. 9: Summary of the "scaled" isovector Dirac and Pauli mean-squared radii from our and RBC/UKQCD Nf = 2 + 1 nucleon 
electromagnetic form- factor calculations [ITI] 



Figure [To] shows the proton Sachs, Dirac and Pauli form factors at the physical pion mass from both quenched and 
dynamical ensembles. The dashed lines are the experimental parametrization as mentioned in earlier subsections. 
We find a significant difference in the results using quenched or dynamical ensembles, most strongly reflected in the 
electric and Dirac form factors. These form factors are better constrained due to the fixed data point at Q^ = 0, 
and the dynamical extrapolations are only a couple of sigma away from the experimental line; the quenched results 
are consistently higher. The magnetic and Pauli form factor extrapolations, on the contrary, fail to reproduce the 
(anomalous) magnetic moments, similar to what we discuss in the magnetic moment subsection, and the extrapolations 
diverge around zero transfer momentum. To our surprise the dynamical and quenched results are consistent for the 
proton Pauli form factor. An examination of the original lattice data to make sure that this is not caused by the 
extrapolation verifies that the similarity comes from the data. This indicates that the Pauli form factor (for the 
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proton) is not very sensitive to whether the QCD vacuum contains fermion degrees of freedom. 

Figure [11] shows a similar plot for the neutron form factors. These include the neutron electric and Dirac form 
factors, although we again caution that there could be significant disconnected contributions that are ignored in this 
calculation. Therefore, it is not fair to compare just the connected contribution with the experimental data; we 
provide it as a guide line to remark on the differences. The smallness of the quenched extrapolation error band is 
due to the large statistics and less-correlated data points, as seen in Fig. |4l We note that as we anticipated, the 
connected-only contribution agrees fairly well at large transfer momentum; this suggests an interpretation assigning 
the remaining discrepancy at low transfer momentum to the omitted disconnected terms. The neutron magnetic and 
Pauli form factors suffer from similar problems to those from the proton, that the fits are less constrained and the 
quenched magnetic form factor has obviously divergent behavior. Also similar to the proton case, the neutron Pauli 
form factor shows indifference to the fermionic vacuum. 

Figure [T^ shows the individual quark and isovector contributions to the Dirac and Pauli form factors. The dashed 
lines are constructed from the proton and neutron electric and magnetic experimental parameterizations form factors, 
as shown in Eq. [151 On the lattice, we can explicitly pick out the individual quark contributions by varying the input 
quark currents. Once again, we see the difference between quenched and dynamical is significant for the Dirac and 
negligible for the Pauli form factors. The dynamical results for both up and down quarks are below the experimental 
reconstruction of the quark contributions, which may vary for future precision neutron form factors. However, the 
linear combinations of up and down for the proton and isovector form factors seem to cancel the difference and become 
agreeable with experiment. The Pauli form factors are within two sigma of the experimental values except in the 
low-Q^ region. 




2' (GeV^) Q" (GeV^) 

FIG. 10: Proton form factors on both quenched (purple) and dynamical (orange) ensembles. The dashed line represents the 
experimental parametrization and the lines with bands are from the lattice calculations. 



D. Transverse Densities 



Following the discussions of Ref. 
in Ref. {491 and deuteron in Ref. H 




26| (and subsequent studies for other hadrons: pion density in Ref. [45|, N-Pu 
see Ref. [43] for a detailed review), we now attempt to translate our extracted 



14 



67 -0.10- 




2 3 

e^ (GeV^) q2 (GeV') 

FIG. 11: Neutron form factors on both quenched (purple) and dynamical (orange) ensembles. The dashed line represents the 
experimental parametrization and the lines with bands are from the lattice calculations. 



nucleon form factors into a description of tlie spatial structure of the nucleon. Although the form factors are typically 
thought of as being the Fourier transforms of the nucleon wavefunction, there is a complication due to the imparted 
momentum at the vertex. Since the incoming and outgoing states have different momentum, their wavefunctions are 
not identical. We can avoid this difhculty by only attempting to describe the spatial structure in the plane transverse 
to an infinite momentum boost. Thus the transverse charge density is defined as the Fourier transform of the form 
factor in such a plane: 



pw-I^^M^y". 



where bold vectors b and q lie in the transverse plane. Equivalently, 

'QdQ 



p{h) = 



27r 



-JoibQ)Fi{Q'), 



(23) 



(24) 



for scalar b, where Jo is a Bessel function. We can perform this integral numerically, using the Fi{Q'^) obtained by 
extrapolating our fit form to the physical pion mass. 

The nature of the Bessel integral, oscillatory and exponentially declining, means that the central core of the 
distribution is most strongly impacted by large-Q^ form factors. If we wish to well characterize this part of the 
nucleon density functions, wc need precise information about form factors in the upper range of transfer momentum. 
We demonstrate this effect by restricting the dynamical data set to the region Q < 2 GcV^ (which is the upper limit 
of transfer momentum for many lattice QCD calculations). Then we compare this result to the density obtained from 
using all the available lattice Q^. Figure [T51 shows in red the effect of omitting the highest Q^ data (above 2 GeV^), 
compared to the blue which uses all Q^ from the dynamical ensembles. The impact is significant in the central core, 
as we anticipated; omitting information about large transfer momenta results in a deviation in the density around 
25%. 

Our results for the transverse charge densities of the proton, neutron, up and down quark contributions, are shown 
in Fig. 1141 along with the same quantity using a parametrization of experimental data from Ref. [l|, |43| . Overall, we 
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FIG. 12: Quark contributions to the form factors on both quenched (purple) and dynamical (orange) ensembles. The dashed 
line represents the experimental parametrization and the lines with bands are from the lattice calculations. 



observe that the dynamical content of the vacuum appears to have a strong influence on the transverse charge densities 
in all cases; we sec large differences between the quenched and dynamical results. In particular, the effect appears 
to increase the density of the core region for quenched ensembles. This difference could be plausibly assigned to the 
presence of a mesonic cloud in the dynamical case, which spreads out the charge distribution. Since the quenched 
ensembles do not include the effects of sea quarks, a nuclcon in that environment docs not couple to baryon-meson 
systems in the usual way. 

In the proton channel, there is quite good agreement between our dynamical and experimental results. This might 
seem suprising, given the omission of the disconnected diagram from our calculation. However, it is a relatively 
small contribution to the proton Dirac form factor; we expect it to be at the level of a percent in this case. Since the 
quenched F\ form factor lies above the experimental one throughout the entire Q^ region we calculated, it accumulates 
a large contribution to the density near the core region. The density becomes smaller than the experimental value at 
larger distances due to the oscillatory nature of the Bessel function. Similar behavior also occurs in other channels. 

The neutron density is expected to have the largest systematics due to the neglect of the disconnected contributions 
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FIG. 13: The transverse charge density p(6) of the proton in the central region with diflerent cuts applied to the -Fi data prior 
to fitting. The lower blue band uses all available Q^ from the 2+1-flavor ensembles; the upper red band uses only data with 



Q < 2 GeV . The experimental result derived from the parametrization of Ref. 



uM 



is shown as a dashed line. 



in this work. Even so, we observe only about a 25% difference between the dynamical results and experiment, 
surprising considering that the disconnected value is on the same order as the neutron form factor itself. Once again, 
the difference becomes much more significant in the quenched case, and the even larger negative core charge density 
is consistent with suppression of mesonic contributions. 

The experimental up and down densities are calculated using linear combinations of the experimental paramctriza- 
tions of the proton and neutron Dirac form factors, as shown in Eqs. [151 The lattice ones are calculated directly 
using inserted vector currents with specific quark flavors. As in the proton case, the disconnected diagrams should 
be minimal. The agreement between the dynamical results and experiment for individual quark contributions is not 
as good as in the case of the proton; in particular, the down-quark density is smaller. It is also notable here that the 
up-quark distribution is more sharply peaked in the core than the down-quark distribution. 

Another kind of transverse density can be derived for the magnetic moment of the nucleons, the transverse mag- 
netization density. The anomalous magnetic moment density is defined as a Fourier transformation of the F2 form 
factor: 



PM(b) 






F2(q2)e 



— iqb 



(25) 



However, Ref. [26] finds it more convenient to consider this in terms of a magnetization density. Therefore, we take a 
derivative with respect to the by direction orthogonal to the magnetic field. This yields 



PAiib) = -by 






= b sin 



'Q^dQ 

2n 



Ji{bQ)F2iQ^), 



(26) 
(27) 



where (p is the angle between b and the magnetic field. Again, we calculate this quantity for the proton, neutron, up- 
and down-quark contributions to the proton and the isovector. We show one-dimensional cuts of the magnetization 
density along the cj) = '"'/S axis of the transverse plane in Fig. [151 We also show the full dependence of the proton and 
neutron magnetization densities across the transverse plane in Fig. 1161 

For the magnetization, we see poorer agreement with experiment than in the case of transverse charge density, 
but better agreement between the quenched and dynamical cases. The Pauli form factors from the dynamical and 
quenched ensembles are relatively similar (unlike the Dirac form factors), and this results in similar magnetization 
densities. The discrepancy in the up-quark channel is strongly reflected in the up-quark dominated proton, despite 
the down-quark contribution having relatively closer agreement with experiment; the converse effect is seen in the 
neutron. 

In Subsec. IIIIBi we found most anomalous magnetic moments from the lattice calculation were roughly 2/3 the 
experimental values; this leads to a discrepancy in the overall scales of the densities. Also the fit forms for the 
Dirac form factors are better constrained in the low transfer-momentum region by conservation of charge, resulting in 
smaller discrepancies at large distances than the magnetization density. To get a better picture of the magnetization 
density, one would also need a better knowledge of the low transfer momentum form factors. Techniques such as 
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FIG. 14: The transverse charge density p{b) of the proton (upper left), neutron (upper right), up (lower left) and down 
(lower right) quark contributions, obtained by integrating the fitted -Fi form at the physical pion mass. Results from the 
dynamical 2+1-flavor ensembles are shown as orange bands. Results from quenched ensembles are shown as purple bands. The 
experimental result derived from the parametrization of Ref [l|, |43] is shown as a dashed line. 

twisted boundary conditions can provide smaller momenta for a particular combination of lattice spacing and box 
size, which should help to ameliorate these discrepancies. 



IV. CONCLUSION 



In this work, we attempted to determine the large-Q^ form factors in lattice QCD. Using operators that have 
various degrees of overlap with boosted states in form-factor calculations significantly improves the signal. Here we 
demonstrated the method using a simple baryon operator with various smearing parameters. Since baryon systems 
typically have worse signal-to-noise ratio than mesons, the method should be applicable to other hadronic systems. 
We analyzed the nucleon ground-state form factor by explicitly including excited-state contributions. This allowed 
us to extract the ground-state signal more precisely, and combined with the use of a longer source-sink separation, 
reduces systematic uncertainty coming from excited-state contamination. 

We checked the most commonly considered observables such as the Dirac and Pauli radii and anomalous magnetic 
moments on our dynamical ensemble, and compared with other Nf = 2-|- 1 calculations, finding reasonable agreement 
within the statistical errors. When we looked at the the dimensionless quantities 4(m'^*)^(r^ 2) in the later section, 
the chiral curvature was significantly reduced in our data; we further found the naive (linear) chiral extrapolation 
recovered the experimental values. Similar behavior was also seen in the RBC/UKQCD lighter pion-mass data points. 
Performing further investigation of such dimensionless quantities with higher signal-to-noise ratios in the lighter pion- 
mass region would be interesting. Nevertheless, the agreement in the low-Q^ behavior among calculations with 
Nf = 2-1-1 demonstrates a nice universality among various fermion actions. 

We extended our study to the high-Q^ region where no lattice-QCD calculations had reached in the past. We found 
significant differences in the Dirac form factor between quenched and dynamical studies, indicating non-negligible 
systematic error due to quenching, while Pauli form factors appear less sensitive to sea-fermion effects. We compared 
our results with an interpolation to experimental results as a function of Q^, although our knowledge of the Q^- 
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FIG. 15: The transverse magnetization density phiib) of the proton (upper left), neutron (upper right), up (lower left) and 
down (lower right) quark contributions, obtained by integrating the fitted F\ form at the physical pion mass. Results from the 
dynamical 2+1-flavor ensembles are shown as orange bands. Results from quenched ensembles are shown as purple bands. The 
experimental result derived from the parametrization of Ref. [l|, |43] is shown as a dashed line. 
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FIG. 16: The two-dimensional dependence of the transverse magnetization density pj\/(b) of the proton (left) and neutron 
(right). Black indicates the lowest values, followed by purple, orange and white in increasing order. For the proton, zero is 
black and values are positive; for the neutron zero is white and values are negative. Figure [15] corresponds to slices through 
the vertical axis. 
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dependence of the G^ form factor remains limited. These experimental lines may change after the collection of future 
precision and large-Q^ neutron form factor data. Our lattice neutron form factors extend as far as 4-6 GcV^ but 
due to the omission of "disconnected" diagrams (which are expected to contribute at O(10~^)), the F" and G^ form 
factors in our calculations suffer from comparable systematic error. However, the other form factors (including those 
for up and down quarks) have relatively large magnitudes, making such O(10~^) systematic uncertainty at the level 
of the statistical error; thus, these should be more reliable compared with experiment. All calculations may be further 
improved with lighter pions and more statistics. 

We look at the transverse charge and magnetization densities using the infinite-momentum frame definition from 
Dirac and Pauli form factors. For these quantities, we see (possibly coincidental) exceptional agreement with exper- 
iment for the charge density of the proton. The level of agreement is particularly striking when compared to the 
quenched results. The neutron charge density agrees less well, as best seen in the isovector charge density, where 
the central density (inside 0.2 fm) greatly exceeds experiment. Since this region is most sensitive to high-Q^ contri- 
butions, large systematics probably exist for both the lattice and experimental measurements. We showed that the 
magnetization densities are much less susceptible to sea-quark effects. In this case, the proton magnetization density 
has greater tension with experiment than the neutron density, but neither is in particularly good agreement. 

We have presented lattice calculations in the large momentum-transfer region up to 4 and 6 GcV . Our currently 
accessible momenta are limited by the available lattice spacing in the calculation, but this will soon be improved as 
finer lattices become available. The signals presented can be further improved by reducing the source-sink separation 
with the same analysis procedure; we plan to proceed in that direction and crosscheck with the calculation done in this 
work. Future generalizations to operators constructed in irreducible representations of the cubic group will probably 
allow us to analyze form factors for radially excited states of the nucleons. 

To reach even higher Q^ regions, we propose performing a numerical step-scaling calculation (for example, Ref. |49|). 
On a small volume with very fine lattice spacing, we can easily reach high momentum. By calculating the step-scaling 
function at overlapping momentum points (or interpolating momentum function) we can reduce the systematic error 
due to finite- volume or lattice discretization artifacts. However, generating several volumes of dynamical Nf = 2 + 1 
lattices requires a large amount of computational resources; we hope such a proposal will become feasible as the 
pctascale computing facilities become available in the near future. 
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